Data-wrangling for spatial dynamics, using "SCTLD_END_exta.csv"
Functions for running and plotting Ripley's K analyses
Figure 7 in paper
library(zoo) # has rollmeans function for calculating AUC
library(igraph)
library(spdep)
library(spatstat)
library(abind)
library(sciplot)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
library(plyr)
library(dplyr)
library(ggpubr)
library(vegan)
library(nlme)
library(car)
library(patchwork)
library(viridis)
getwd()
[1] "/Users/saradellwilliams/Dropbox/My_Publications/SpatEpiSCTLD/SCTLDepizootiology_lowerFLkeys/SpatialAnalyses"
my.data<-read.csv("SCTLD_END_exta.csv") ## this file is shortened and doesnt included the treatment sites
#move to a long format so that every row is now an observation of a colony at a single timepoint
data_long <- gather(my.data, key=timept, value=state, X5.1.18:X12.6.19,factor_key = TRUE)
data_long<-data_long[,-c(8,9)] #drop column 8:10, which is glom and not needed anymore
#head(data_long)
#make the health states factors
data_long$state<-factor(data_long$state,levels=c("Healthy","SCTLD","Dead","Unknown"))
#revalue the timepoints
data_long$timept<-revalue(data_long$timept, c( "X5.1.18"="05-10-18","X6.1.18"="06-01-18","X6.21.18"="06-21-18","X7.16.18"="07-16-18","X8.17.18"="08-17-18","X10.30.18"="10-30-18", "X11.9.18"="11-09-18", "X11.29.18"="11-29-18","X12.13.18"="12-13-18","X1.4.19"="01-04-19","X1.18.19"="01-18-19","X2.8.19"="02-08-19","X3.4.19"="03-04-19","X3.21.19"="03-21-19","X4.11.19"="04-11-19","X5.2.19"="05-02-19","X5.16.19"="05-16-19","X5.28.19"="05-28-19","X6.13.19"="06-13-19","X7.1.19"="07-01-19","X7.22.19"="07-22-19","X8.16.19"="08-16-19","X9.17.19"="09-17-19","X10.14.19"="10-14-19","X11.12.19"="11-12-19","X12.6.19"="12-06-19"))
head(data_long)
Site Plot Sps Max_width Coral_ID coords_x coords_y days_dis
1 1 23 DLAB 11 1_p23_t1_s0_c1_DLAB 0.448071 0.2361312 0
2 1 23 SINT 22 1_p23_t1_s0_c2_SINT 0.789318 0.4790354 0
3 1 23 SSID 25 1_p23_t1_s0_c3_SSID 0.833828 0.6278998 0
4 1 23 PPOR 11 1_p23_t1_s0_c4_PPOR 0.225519 1.9267537 0
5 1 23 DSTO 18 1_p23_t1_s0_c5_DSTO 0.462908 2.2250339 164
6 1 23 CNAT 11 1_p23_t1_s0_c6_CNAT 0.151335 2.2052883 0
weeks_dis timept state
1 0.00000 05-10-18 Healthy
2 0.00000 05-10-18 Healthy
3 0.00000 05-10-18 Healthy
4 0.00000 05-10-18 Healthy
5 23.42857 05-10-18 Healthy
6 0.00000 05-10-18 Healthy
summary(data_long$Plot)
Min. 1st Qu. Median Mean 3rd Qu. Max.
23.00 25.00 28.00 33.57 45.00 47.00
##keep only susceptible species
data_filtered<-data_long%>%
filter(Sps!="AAGA",Sps!="ACER",Sps!="ATEN",Sps!="EFAS",Sps!="MANG",Sps!="MMEA",Sps!="MYCE",Sps!="OCUL",Sps!="ODIF",Sps!="PAST",Sps!="PDIV",Sps!="PPOR",Sps!="SRAD")%>%
droplevels()
Plot45 <-data_filtered %>%
dplyr::filter(coords_x!="NA",Plot==45)
Plot47 <- data_filtered %>%
dplyr::filter(coords_x!="NA",Plot==47)
Into timepoints of disease colonies per plot
p47<-as.data.frame(Plot47)
p45<-as.data.frame(Plot45)
#str(p47)
#now just disease
p47dis<-as.data.frame(filter(Plot47,state=="SCTLD"))
p45dis<-as.data.frame(filter(Plot45,state=="SCTLD"))
#split 47 disease tps
#start on March 4 because that's when disease cases >5 for each
p47dis_tp1<-(subset(p47dis,timept=="03-04-19",select=c(coords_x,coords_y)))
p47dis_tp2<-(subset(p47dis,timept=="03-21-19",select=c(coords_x,coords_y)))
p47dis_tp3<-(subset(p47dis,timept=="04-11-19",select=c(coords_x,coords_y)))
p47dis_tp4<-(subset(p47dis,timept=="05-02-19",select=c(coords_x,coords_y)))
p47dis_tp5<-(subset(p47dis,timept=="05-16-19",select=c(coords_x,coords_y)))
p47dis_tp6<-(subset(p47dis,timept=="05-28-19",select=c(coords_x,coords_y)))
p47dis_tp7<-(subset(p47dis,timept=="06-13-19",select=c(coords_x,coords_y)))
p47dis_tp8<-(subset(p47dis,timept=="07-01-19",select=c(coords_x,coords_y)))
p47dis_tp9<-(subset(p47dis,timept=="07-22-19",select=c(coords_x,coords_y)))
p47dis_tp10<-(subset(p47dis,timept=="08-16-19",select=c(coords_x,coords_y)))
p47dis_tp11<-(subset(p47dis,timept=="09-17-19",select=c(coords_x,coords_y)))
p47dis_tp12<-(subset(p47dis,timept=="10-14-19",select=c(coords_x,coords_y)))
p47dis_tp13<-(subset(p47dis,timept=="11-12-19",select=c(coords_x,coords_y)))
p47dis_tp14<-(subset(p47dis,timept=="12-06-19",select=c(coords_x,coords_y)))
#split 45 disease tps
p45dis_tp1<-(subset(p45dis,timept=="03-04-19",select=c(coords_x,coords_y)))
p45dis_tp2<-(subset(p45dis,timept=="03-21-19",select=c(coords_x,coords_y)))
p45dis_tp3<-(subset(p45dis,timept=="04-11-19",select=c(coords_x,coords_y)))
p45dis_tp4<-(subset(p45dis,timept=="05-02-19",select=c(coords_x,coords_y)))
p45dis_tp5<-(subset(p45dis,timept=="05-16-19",select=c(coords_x,coords_y)))
p45dis_tp6<-(subset(p45dis,timept=="05-28-19",select=c(coords_x,coords_y)))
p45dis_tp7<-(subset(p45dis,timept=="06-13-19",select=c(coords_x,coords_y)))
p45dis_tp8<-(subset(p45dis,timept=="07-01-19",select=c(coords_x,coords_y)))
p45dis_tp9<-(subset(p45dis,timept=="07-22-19",select=c(coords_x,coords_y)))
p45dis_tp10<-(subset(p45dis,timept=="08-16-19",select=c(coords_x,coords_y)))
p45dis_tp11<-(subset(p45dis,timept=="09-17-19",select=c(coords_x,coords_y)))
p45dis_tp12<-(subset(p45dis,timept=="10-14-19",select=c(coords_x,coords_y)))
p45dis_tp13<-(subset(p45dis,timept=="11-12-19",select=c(coords_x,coords_y)))
p45dis_tp14<-(subset(p45dis,timept=="12-06-19",select=c(coords_x,coords_y)))
#get all cols at one timepoint
p47.all<-(subset(p47,timept=="06-13-19",select=c(coords_x,coords_y)))
#p47.all
p45.all<-(subset(p45,timept=="06-13-19",select=c(coords_x,coords_y)))
#Runs and Plots results of Ripleys K
myfunc.ppp.Kest<-function(data,w,name){
data.ppp <- as.ppp(data,w)
# plot(data.ppp)
lhatK <-Kest(data.ppp)
#dividing by pi and taking the square root, linearizes the Ripley's value
plot(lhatK$r,sqrt(lhatK$iso/pi)-lhatK$r,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",col="green",ylab="L(r)",lty=1,main=name)
# plot(lhatK$r,lhatK$iso,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",col="green",ylab="L(r)",lty=1,main=name)
return(c(data.ppp,lhatK)) #return the lhatK and ppp for use later if needed
}
#### Function to make null distribution Plot
#allpts & dispts should be my normal dataframes
myfunc.nulldistrKest<-function(allpts,dispts,w,n,name){
#get info for all
all<-as.data.frame(t(allpts))
allX<-t(all)
all.spp <- unique(as.ppp(allX,w))
#print(anyDuplicated(all.spp))
all.lhat <-Kest(all.spp)
#get info for dis
dis<-as.data.frame(t(dispts))
disX<-t(dis)
dis.spp <- unique(as.ppp(disX,w))
#print(anyDuplicated(dis.spp))
dis.lhat <-Kest(dis.spp)
plot(dis.lhat$r,sqrt(dis.lhat$iso/pi)-dis.lhat$r,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",ylab="L(r)",lty=0,pch=NA,main=name)
for (i in 1:1000){
all.samp<-sample(all,n,replace=FALSE)# n=total number of disease points to sample from all point in the plot
x <-t(all.samp)
all.samp.spp <- unique(as.ppp(x,w))
#print(anyDuplicated(all.samp.spp))
all.samp.lhat <-Kest(all.samp.spp)
lines(all.samp.lhat$r,sqrt(all.samp.lhat$iso/pi)-all.samp.lhat$r,xlab="r (meters)",col="gray",lty=1,lwd=1)
}
lines(all.lhat$r,sqrt(all.lhat$theo/pi)-all.lhat$r,lwd=1,lty=8,col="black")
lines(dis.lhat$r,sqrt(dis.lhat$iso/pi)-dis.lhat$r,lwd=1,lty=1,col="red")
Hmisc::minor.tick(nx=10, ny=5, tick.ratio=0.5)
}
#### Functions to get difference in L(R) between actual and null distributions ######
### Function to get top of null distribution #####
top_of_RKnull<-function(allpts,w,n){
all<-as.data.frame(t(allpts))
allX<-t(all)
all.spp <- unique(as.ppp(allX,w))
all.rk <-Kest(all.spp)
all.lhat_adj<-sqrt(all.rk$iso/pi)-all.rk$r
plot(all.rk$r,sqrt(all.rk$iso/pi)-all.rk$r,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",ylab="L(r)",lty=1)
for (i in 1:1000){
all.samp<-sample(all,n,replace=FALSE)# n=total number of disease points to sample from all point in the plot
x <-t(all.samp)
all.samp.spp <- unique(as.ppp(x,w))
#print(anyDuplicated(all.samp.spp))
all.samp.rk <-Kest(all.samp.spp)
l_adj<-sqrt(all.samp.rk$iso/pi)-all.samp.rk$r
all.lhat_adj<-cbind(all.lhat_adj,l_adj)
lines(all.samp.rk$r,sqrt(all.samp.rk$iso/pi)-all.samp.rk$r,xlab="r (meters)",col="gray",lty=1,lwd=1)
}
##radius stays same each time
linemaxs<-colSums(all.lhat_adj)
maxline<-max(linemaxs)
best_index<-which(linemaxs==max(colSums(all.lhat_adj)))
null_line<-all.lhat_adj[,best_index]
return(null_line)
}
### Function to get difference between top of null and disease line #####
lhat_diffs<-function(null_line,dispts,w,n){
dis<-as.data.frame(t(dispts))
disX<-t(dis)
dis.spp <- unique(as.ppp(disX,w))
dis.rk <-Kest(dis.spp)
dis.lhat_adj<-sqrt(dis.rk$iso/pi)-dis.rk$r
lhat_dif<-dis.lhat_adj-null_line
plot(dis.rk$r,lhat_dif,type="l")
abline(h=0)
results<-cbind(dis.rk$r,dis.lhat_adj,null_line,lhat_dif)
return(results)
}
### Function to get top of null distribution WITH PLOT #####
top_of_RKnull_keepforplot<-function(allpts,w,n){
all<-as.data.frame(t(allpts))
allX<-t(all)
all.spp <- unique(as.ppp(allX,w))
all.rk <-Kest(all.spp)
all.lhat_adj<-sqrt(all.rk$iso/pi)-all.rk$r
my.mat<-matrix(nrow=513,ncol=1000)
plot(all.rk$r,sqrt(all.rk$iso/pi)-all.rk$r,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",ylab="L(r)",lty=1)
for (i in 1:1000){
all.samp<-sample(all,n,replace=FALSE)# n=total number of disease points to sample from all point in the plot
x <-t(all.samp)
all.samp.spp <- unique(as.ppp(x,w))
#print(anyDuplicated(all.samp.spp))
all.samp.rk <-Kest(all.samp.spp)
l_adj<-sqrt(all.samp.rk$iso/pi)-all.samp.rk$r
all.lhat_adj<-cbind(all.lhat_adj,l_adj)
lines(all.samp.rk$r,sqrt(all.samp.rk$iso/pi)-all.samp.rk$r,xlab="r (meters)",col="gray",lty=1,lwd=1)
my.mat[,i]<-sqrt(all.samp.rk$iso/pi)-all.samp.rk$r
}
##radius stays same each time
linemaxs<-colSums(all.lhat_adj)
maxline<-max(linemaxs)
best_index<-which(linemaxs==max(colSums(all.lhat_adj)))
null_line<-all.lhat_adj[,best_index]
my.df<-data.frame(my.mat)
results.df<-cbind(null_line,my.mat)
return(results.df)
}
##### func for using prev func output for plot #######
myfunc.nullplot2<-function(null_line,dispts,allpts,w,n,name){
#get info for all
all<-as.data.frame(t(allpts))
allX<-t(all)
all.spp <- unique(as.ppp(allX,w))
#print(anyDuplicated(all.spp))
all.lhat <-Kest(all.spp)
#get info for dis
dis<-as.data.frame(t(dispts))
disX<-t(dis)
dis.spp <- unique(as.ppp(disX,w))
#print(anyDuplicated(dis.spp))
dis.lhat <-Kest(dis.spp)
plot(dis.lhat$r,sqrt(dis.lhat$iso/pi)-dis.lhat$r,cex.axis=1.5,cex.lab=1.5,xlab="r (m)",ylab="L(r)",lty=0,pch=NA,main=name)
for (i in 1:1000){
y<-null_line[,i+1]
lines(all.lhat$r[0:nrow(null_line)],y,xlab="r (meters)",col="gray",lty=1,lwd=1)
}
lines(all.lhat$r,sqrt(all.lhat$theo/pi)-all.lhat$r,lwd=1,lty=8,col="black")
lines(dis.lhat$r,sqrt(dis.lhat$iso/pi)-dis.lhat$r,lwd=1,lty=1,col="red")
Hmisc::minor.tick(nx=10, ny=5, tick.ratio=0.5)
}
w <- owin(c(-20,20),c(-20,20))
n <- 100
allpts<-p47.all
dispts<-p47dis_tp6
name<-"Quadrat 47 5/28/19"
myfunc.nulldistrKest(allpts,dispts,w,n,name)
Since the null distribution is random each time, I saved the null distribution for each plot once and used the same null distribution for each time point
p47.nl1<-read.csv("plot47RKnulldist.csv",colClasses = "numeric")
myfunc.nullplot2(p47.nl1,p47dis_tp6,p47.all ,w,n=100,"Plot 47, May 28th 2019")
p45.nl1<-read.csv("plot45RKnulldist.csv",colClasses = "numeric")
myfunc.nullplot2(p45.nl1,p45dis_tp6,p45.all ,w,n=100,"Plot 45, May 28th 2019")
data contain duplicated points
#Saved just the top of the null distribution for each each quadrat which is used for getting values.
p47.nl1<-read.csv("topofthenull_plot47.csv",header=T)
p45.nl1<-read.csv("topofthenull_plot45.csv",header=T)
w <- owin(c(-20,20),c(-20,20))
n <- 100
head(p45.nl1)
X x
1 1 0.00000000
2 2 -0.01953125
3 3 -0.03906250
4 4 -0.05859375
5 5 -0.07812500
6 6 -0.09765625
w <- owin(c(-20,20),c(-20,20))
n <- 100
# GET DIFS BETWEEN ALL TP DIS AND TOP OF NULL
p45ld.tp1<-lhat_diffs(p45.nl1[,2],p45dis_tp1,w,n)
p45ld.tp2<-lhat_diffs(p45.nl1[,2],p45dis_tp2,w,n)
p45ld.tp3<-lhat_diffs(p45.nl1[,2],p45dis_tp3,w,n)
p45ld.tp4<-lhat_diffs(p45.nl1[,2],p45dis_tp4,w,n)
p45ld.tp5<-lhat_diffs(p45.nl1[,2],p45dis_tp5,w,n)
p45ld.tp6<-lhat_diffs(p45.nl1[,2],p45dis_tp6,w,n)
p45ld.tp7<-lhat_diffs(p45.nl1[,2],p45dis_tp7,w,n)
p45ld.tp8<-lhat_diffs(p45.nl1[,2],p45dis_tp8,w,n)
p45ld.tp9<-lhat_diffs(p45.nl1[,2],p45dis_tp9,w,n)
p45ld.tp10<-lhat_diffs(p45.nl1[,2],p45dis_tp10,w,n)
p45ld.tp11<-lhat_diffs(p45.nl1[,2],p45dis_tp11,w,n)
p45ld.tp12<-lhat_diffs(p45.nl1[,2],p45dis_tp12,w,n)
p45ld.tp13<-lhat_diffs(p45.nl1[,2],p45dis_tp13,w,n)
p45ld.tp14<-lhat_diffs(p45.nl1[,2],p45dis_tp14,w,n)
### now for plot 47
p47ld.tp1<-lhat_diffs(p47.nl1[,2],p47dis_tp1,w,n)
p47ld.tp2<-lhat_diffs(p47.nl1[,2],p47dis_tp2,w,n)
p47ld.tp3<-lhat_diffs(p47.nl1[,2],p47dis_tp3,w,n)
p47ld.tp4<-lhat_diffs(p47.nl1[,2],p47dis_tp4,w,n)
p47ld.tp5<-lhat_diffs(p47.nl1[,2],p47dis_tp5,w,n)
p47ld.tp6<-lhat_diffs(p47.nl1[,2],p47dis_tp6,w,n)
p47ld.tp7<-lhat_diffs(p47.nl1[,2],p47dis_tp7,w,n)
p47ld.tp8<-lhat_diffs(p47.nl1[,2],p47dis_tp8,w,n)
p47ld.tp9<-lhat_diffs(p47.nl1[,2],p47dis_tp9,w,n)
p47ld.tp10<-lhat_diffs(p47.nl1[,2],p47dis_tp10,w,n)
p47ld.tp11<-lhat_diffs(p47.nl1[,2],p47dis_tp11,w,n)
p47ld.tp12<-lhat_diffs(p47.nl1[,2],p47dis_tp12,w,n)
p47ld.tp13<-lhat_diffs(p47.nl1[,2],p47dis_tp13,w,n)
p47ld.tp14<-lhat_diffs(p47.nl1[,2],p47dis_tp14,w,n)
p47.ldiffs<-(cbind(p47ld.tp1[,1],p47ld.tp1[,4],p47ld.tp2[,4],p47ld.tp3[,4],p47ld.tp4[,4],p47ld.tp5[,4],p47ld.tp6[,4],p47ld.tp7[,4],p47ld.tp8[,4],p47ld.tp9[,4],p47ld.tp10[,4],p47ld.tp11[,4],p47ld.tp12[,4],p47ld.tp13[,4],p47ld.tp14[,4]))
p47.ldiffs<-data.frame(p47.ldiffs)
colnames(p47.ldiffs)<-c("r","03-04-19","03-21-19","04-11-19","05-02-19","05-16-19","05-28-19","06-13-19","07-01-19","07-22-19","08-16-19","09-17-19","10-14-19","11-12-19","12-06-19")
#now plot 45
p45.ldiffs<-(cbind(p45ld.tp1[,1],p45ld.tp1[,4],p45ld.tp2[,4],p45ld.tp3[,4],p45ld.tp4[,4],p45ld.tp5[,4],p45ld.tp6[,4],p45ld.tp7[,4],p45ld.tp8[,4],p45ld.tp9[,4],p45ld.tp10[,4],p45ld.tp11[,4],p45ld.tp12[,4],p45ld.tp13[,4],p45ld.tp14[,4]))
p45.ldiffs<-data.frame(p45.ldiffs)
colnames(p45.ldiffs)<-c("r","03-04-19","03-21-19","04-11-19","05-02-19","05-16-19","05-28-19","06-13-19","07-01-19","07-22-19","08-16-19","09-17-19","10-14-19","11-12-19","12-06-19")
#nrow(p45.ldiffs)
For Quadrat/Plot 45, there are two peak clustering radii
#look over the 25 point range around the max peak to get the best value of it.
maxdiff.low<-matrix(nrow=25,ncol=14)
radius.low<-matrix(nrow=25,ncol=14)
maxdiff.high<-matrix(nrow=25,ncol=14)
radius.high<-matrix(nrow=25,ncol=14)
###lower clusters
p45.ldiffs.copy<-p45.ldiffs[1:205,]
for (i in 2:ncol(p45.ldiffs.copy)){
p45.ldiffs.copy<-p45.ldiffs[1:205,]
for (j in 1:25){
maxdiff.low[j,i-1]<-max(p45.ldiffs.copy[,i])
index<-(which(p45.ldiffs.copy[,i]==max(p45.ldiffs.copy[,i])))
#print(index[1])
crad<-p45.ldiffs.copy$r[index[1]]
#print(crad)
radius.low[j,i-1]<-crad
p45.ldiffs.copy<-p45.ldiffs.copy[-index,]
}
}
p45.ldiffs.copy<-p45.ldiffs[205:513,]
### upper clusters
for (i in 2:ncol(p45.ldiffs.copy)){
p45.ldiffs.copy<-p45.ldiffs[205:513,]
for (j in 1:25){
maxdiff.high[j,i-1]<-max(p45.ldiffs.copy[,i])
index<-(which(p45.ldiffs.copy[,i]==max(p45.ldiffs.copy[,i])))
#print(index[1])
crad<-p45.ldiffs.copy$r[index[1]]
#print(crad)
radius.high[j,i-1]<-crad
p45.ldiffs.copy<-p45.ldiffs.copy[-index,]
}
}
radmat.p45.low<-radius.low
maxdiffmat.p45.low<-maxdiff.low
radmat.p45.high<-radius.high
maxdiffmat.p45.high<-maxdiff.high
maxdiffmat.p45.low[,3]
[1] 0.0000000 -0.1710455 -0.2215138 -0.2491504 -0.2779359 -0.2811128 -0.3023020
[8] -0.3451862 -0.3556473 -0.3566962 -0.3837333 -0.4020600 -0.4224172 -0.4396258
[15] -0.4536255 -0.4575955 -0.4627132 -0.4710339 -0.4716195 -0.4869890 -0.5331168
[22] -0.5505893 -0.5520653 -0.5555755 -0.5561665
colMeans(radmat.p45.high)
[1] 4.316406 4.348438 7.422656 5.320312 8.236719 4.397656 6.658594 5.912500
[9] 5.912500 6.855469 6.592969 6.625781 5.882812 7.741406
colMeans(radmat.p45.low)
[1] 1.0625000 3.7296875 1.9507812 2.3085938 1.2609375 0.6078125 1.5859375
[8] 1.7031250 1.7031250 1.8945312 1.8945312 1.8210937 3.0078125 1.7875000
For quadrat/plot 47, there is just one peak clustering radius
maxdiff<-matrix(nrow=50,ncol=14)
radius<-matrix(nrow=50,ncol=14)
p47.ldiffs.copy<-p47.ldiffs
#p47.ldiffs.copy
p47.ldiffs.copy<-p47.ldiffs
for (i in 2:ncol(p47.ldiffs.copy)){
p47.ldiffs.copy<-p47.ldiffs
for (j in 1:50){
maxdiff[j,i-1]<-max(p47.ldiffs.copy[,i])
index<-(which(p47.ldiffs.copy[,i]==max(p47.ldiffs.copy[,i])))
#print(index[1])
crad<-p47.ldiffs.copy$r[index[1]]
#print(crad)
radius[j,i-1]<-crad
p47.ldiffs.copy<-p47.ldiffs.copy[-index,]
}
}
radmat.p47<-radius
maxdiffmat.p47<-maxdiff
colMeans(radmat.p47)
[1] 1.5277344 1.7343750 3.8015625 4.5398438 3.2835938 3.6539062 3.1726563
[8] 3.0089844 3.0089844 2.6308594 2.5714844 2.5925781 2.5753906 0.9734375
upR45<-colMeans(radmat.p45.high)
lowR45<-colMeans(radmat.p45.low)
R47<-colMeans(radmat.p47)
datepts<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
plot(datepts,upR45,ylim=c(0,10),t="l",col="green")
lines(datepts,lowR45,col="blue")
lines(datepts,R47,col="black",pch=2)
R.results<-cbind(c("03-04-19","03-21-19","04-11-19","05-02-19","05-16-19","05-28-19","06-13-19","07-01-19","07-22-19","08-16-19","09-17-19","10-14-19","11-12-19","12-06-19"),
upR45,lowR45,R47)
plot 47 first
top_range<-c()
bottom_range<-c()
#so we want to get the range of significant R
for (j in 2:ncol(p47.ldiffs)){
#print(j)
#print(R47[j-1])
#the top half
top_start<-which(p47.ldiffs[,1]>=R47[j-1])[1]
#print(top_start)
top_ranges<-c()
for(i in top_start:nrow(p47.ldiffs)-1){
if(p47.ldiffs[i,j]>=0 & p47.ldiffs[i+1,j]<=0){
#print ("intercept")
top_ranges<-c(top_ranges,p47.ldiffs[i,1])
}
}
if(is.null(top_ranges)){
#print("inf")
top_ranges<-10
}
top_range<-c(top_range,top_ranges[1])
#print(top_range)
#the bottom half
bottom_start<-which(p47.ldiffs[,1]<=R47[j-1])
bottom_start<-bottom_start[length(bottom_start)]
bottom_start
bottom_ranges<-c()
for(i in 1:bottom_start){
if(p47.ldiffs[i,j]<=0 & p47.ldiffs[i+1,j]>=0){
#print ("intercept")
bottom_ranges<-c(bottom_ranges,p47.ldiffs[i,1])
}
}
bottom_range<-c(bottom_range, bottom_ranges[length(bottom_ranges)])
#print(bottom_range)
}
length(bottom_range)
[1] 14
length(R47)
[1] 14
length(top_range)
[1] 14
top_range
[1] 2.929688 2.753906 6.132812 10.000000 10.000000 9.375000 8.906250
[8] 9.042969 9.042969 8.378906 8.476562 9.492188 8.554688 1.542969
R47
[1] 1.5277344 1.7343750 3.8015625 4.5398438 3.2835938 3.6539062 3.1726563
[8] 3.0089844 3.0089844 2.6308594 2.5714844 2.5925781 2.5753906 0.9734375
R47_withrange<-cbind(bottom_range,R47,top_range)
now plot 45, higher R
top_range<-c()
bottom_range<-c()
#so we want to get the range of significant R
for (j in 2:ncol(p45.ldiffs)){
#print(j)
#print(R47[j-1])
#the top half
top_start<-which(p45.ldiffs[,1]>=upR45[j-1])[1]
#print(top_start)
top_ranges<-c()
for(i in top_start:nrow(p45.ldiffs)-1){
if(p45.ldiffs[i,j]>=0 & p45.ldiffs[i+1,j]<=0){
#print ("intercept")
top_ranges<-c(top_ranges,p45.ldiffs[i,1])
}
}
if(is.null(top_ranges)){
#print("inf")
top_ranges<-10
}
top_range<-c(top_range,top_ranges[1])
#print(top_range)
#the bottom half
bottom_start<-which(p45.ldiffs[,1]<=upR45[j-1])
bottom_start<-bottom_start[length(bottom_start)]
bottom_start
bottom_ranges<-c()
for(i in 1:bottom_start){
if(p45.ldiffs[i,j]<=0 & p45.ldiffs[i+1,j]>=0){
#print ("intercept")
bottom_ranges<-c(bottom_ranges,p45.ldiffs[i,1])
}
}
bottom_range<-c(bottom_range, bottom_ranges[length(bottom_ranges)])
#print(bottom_range)
}
length(bottom_range)
[1] 14
length(upR45)
[1] 14
length(top_range)
[1] 14
upR45_withrange<-cbind(bottom_range,upR45,top_range)
upR45_withrange
bottom_range upR45 top_range
[1,] 3.6914062 4.316406 6.757812
[2,] 2.4804688 4.348438 10.000000
[3,] 7.3437500 7.422656 7.519531
[4,] 0.7226562 5.320312 9.492188
[5,] 6.1328125 8.236719 10.000000
[6,] 3.6718750 4.397656 4.511719
[7,] 5.0781250 6.658594 10.000000
[8,] 5.1562500 5.912500 10.000000
[9,] 5.1562500 5.912500 10.000000
[10,] 6.1914062 6.855469 10.000000
[11,] 5.3320312 6.592969 10.000000
[12,] 5.2734375 6.625781 10.000000
[13,] 5.3125000 5.882812 10.000000
[14,] 7.1679688 7.741406 9.082031
#row 5 may not be good
Plot 45, lower R
top_range<-c()
bottom_range<-c()
#so we want to get the range of significant R
for (j in 2:ncol(p45.ldiffs)){
#print(j)
#print(R47[j-1])
#the top half
top_start<-which(p45.ldiffs[,1]>=lowR45[j-1])[1]
#print(top_start)
top_ranges<-c()
for(i in top_start:nrow(p45.ldiffs)-1){
if(p45.ldiffs[i,j]>=0 & p45.ldiffs[i+1,j]<=0){
#print ("intercept")
top_ranges<-c(top_ranges,p45.ldiffs[i,1])
}
}
if(is.null(top_ranges)){
#print("inf")
top_ranges<-10
}
top_range<-c(top_range,top_ranges[1])
#print(top_range)
#the bottom half
bottom_start<-which(p45.ldiffs[,1]<=lowR45[j-1])
bottom_start<-bottom_start[length(bottom_start)]
bottom_start
bottom_ranges<-c()
for(i in 1:bottom_start){
if(p45.ldiffs[i,j]<=0 & p45.ldiffs[i+1,j]>=0){
#print ("intercept")
bottom_ranges<-c(bottom_ranges,p45.ldiffs[i,1])
}
}
bottom_range<-c(bottom_range, bottom_ranges[length(bottom_ranges)])
#print(bottom_range)
}
length(bottom_range)
[1] 14
length(lowR45)
[1] 14
length(top_range)
[1] 14
lowR45_withrange<-cbind(bottom_range,lowR45,top_range)
lowR45_withrange
bottom_range lowR45 top_range
[1,] 0.8007812 1.0625000 2.285156
[2,] 2.4804688 3.7296875 10.000000
[3,] 0.1171875 1.9507812 6.816406
[4,] 0.7226562 2.3085938 9.492188
[5,] 0.2539062 1.2609375 2.656250
[6,] 0.2539062 0.6078125 1.113281
[7,] 1.2304688 1.5859375 3.339844
[8,] 1.2304688 1.7031250 2.480469
[9,] 1.2304688 1.7031250 2.480469
[10,] 1.6406250 1.8945312 2.285156
[11,] 1.3085938 1.8945312 2.480469
[12,] 1.6406250 1.8210937 4.980469
[13,] 2.7148438 3.0078125 4.023438
[14,] 0.7226562 1.7875000 3.085938
#row 5 may not be good
p47.ldiffsc<-p47.ldiffs
p47.ldiffsc[p47.ldiffsc<0]<-0
p47.tp1AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,2],2))
p47.tp2AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,3],2))
p47.tp3AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,4],2))
p47.tp4AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,5],2))
p47.tp5AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,6],2))
p47.tp6AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,7],2))
p47.tp7AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,8],2))
p47.tp8AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,9],2))
p47.tp9AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,10],2))
p47.tp10AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,11],2))
p47.tp11AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,12],2))
p47.tp12AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,13],2))
p47.tp13AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,14],2))
p47.tp14AUC <- sum(diff(p47.ldiffsc[,1])*rollmean(p47.ldiffsc[,15],2))
p47.AUC<-rbind(p47.tp1AUC,p47.tp2AUC,p47.tp3AUC,p47.tp4AUC,p47.tp5AUC,p47.tp6AUC,p47.tp7AUC,p47.tp8AUC,p47.tp9AUC,p47.tp10AUC,p47.tp11AUC,p47.tp12AUC,p47.tp13AUC,p47.tp14AUC)
p47.AUC
[,1]
p47.tp1AUC 5.1360509
p47.tp2AUC 1.2491401
p47.tp3AUC 1.1287722
p47.tp4AUC 9.5481176
p47.tp5AUC 7.3403452
p47.tp6AUC 11.0784857
p47.tp7AUC 11.4838203
p47.tp8AUC 9.7677461
p47.tp9AUC 9.7677461
p47.tp10AUC 9.8531471
p47.tp11AUC 11.5345669
p47.tp12AUC 15.1845803
p47.tp13AUC 5.9864486
p47.tp14AUC 0.4679849
p45.ldiffsc<-p45.ldiffs
p45.ldiffsc[p45.ldiffsc<0]<-0
p45.tp1AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,2],2))
p45.tp2AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,3],2))
p45.tp3AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,4],2))
p45.tp4AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,5],2))
p45.tp5AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,6],2))
p45.tp6AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,7],2))
p45.tp7AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,8],2))
p45.tp8AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,9],2))
p45.tp9AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,10],2))
p45.tp10AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,11],2))
p45.tp11AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,12],2))
p45.tp12AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,13],2))
p45.tp13AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,14],2))
p45.tp14AUC <- sum(diff(p45.ldiffsc[,1])*rollmean(p45.ldiffsc[,15],2))
p45.AUC<-rbind(p45.tp1AUC,p45.tp2AUC,p45.tp3AUC,p45.tp4AUC,p45.tp5AUC,p45.tp6AUC,p45.tp7AUC,p45.tp8AUC,p45.tp9AUC,p45.tp10AUC,p45.tp11AUC,p45.tp12AUC,p45.tp13AUC,p45.tp14AUC)
p45.AUC
[,1]
p45.tp1AUC 12.79715782
p45.tp2AUC 16.58548352
p45.tp3AUC 0.01276313
p45.tp4AUC 10.76357877
p45.tp5AUC 4.04014681
p45.tp6AUC 3.07665302
p45.tp7AUC 8.10033025
p45.tp8AUC 7.05207728
p45.tp9AUC 7.05207728
p45.tp10AUC 3.35213051
p45.tp11AUC 5.40229753
p45.tp12AUC 4.02300920
p45.tp13AUC 3.98433812
p45.tp14AUC 0.72366597
dates<-c("03-04-19","03-21-19","04-11-19","05-02-19","05-16-19","05-28-19","06-13-19","07-01-19","07-22-19","08-16-19","09-17-19","10-14-19","11-12-19","12-06-19")
auc.df<-data.frame(cbind(dates,p45.AUC[,1],p47.AUC[,1]))
colnames(auc.df)<-c("dates","p45auc","p47auc")
rownames(auc.df)<-c()
Rk.results<-cbind(auc.df,upR45_withrange,lowR45_withrange,R47_withrange)
Rk.results
dates p45auc p47auc bottom_range upR45 top_range
1 03-04-19 12.7971578186068 5.13605089602507 3.6914062 4.316406 6.757812
2 03-21-19 16.5854835240743 1.24914007438787 2.4804688 4.348438 10.000000
3 04-11-19 0.0127631255303744 1.12877218748147 7.3437500 7.422656 7.519531
4 05-02-19 10.7635787727553 9.5481176180383 0.7226562 5.320312 9.492188
5 05-16-19 4.04014681315698 7.34034521885237 6.1328125 8.236719 10.000000
6 05-28-19 3.07665301853824 11.0784857151001 3.6718750 4.397656 4.511719
7 06-13-19 8.10033024835652 11.4838203028536 5.0781250 6.658594 10.000000
8 07-01-19 7.0520772791742 9.76774607153853 5.1562500 5.912500 10.000000
9 07-22-19 7.0520772791742 9.76774607153853 5.1562500 5.912500 10.000000
10 08-16-19 3.35213050996111 9.8531471166574 6.1914062 6.855469 10.000000
11 09-17-19 5.40229752523318 11.5345668776599 5.3320312 6.592969 10.000000
12 10-14-19 4.02300919883857 15.1845803014556 5.2734375 6.625781 10.000000
13 11-12-19 3.9843381193459 5.9864486039555 5.3125000 5.882812 10.000000
14 12-06-19 0.723665974475861 0.467984852389785 7.1679688 7.741406 9.082031
bottom_range lowR45 top_range bottom_range R47 top_range
1 0.8007812 1.0625000 2.285156 0.80078125 1.5277344 2.929688
2 2.4804688 3.7296875 10.000000 1.30859375 1.7343750 2.753906
3 0.1171875 1.9507812 6.816406 2.44140625 3.8015625 6.132812
4 0.7226562 2.3085938 9.492188 1.01562500 4.5398438 10.000000
5 0.2539062 1.2609375 2.656250 0.33203125 3.2835938 10.000000
6 0.2539062 0.6078125 1.113281 0.60546875 3.6539062 9.375000
7 1.2304688 1.5859375 3.339844 0.33203125 3.1726563 8.906250
8 1.2304688 1.7031250 2.480469 0.33203125 3.0089844 9.042969
9 1.2304688 1.7031250 2.480469 0.33203125 3.0089844 9.042969
10 1.6406250 1.8945312 2.285156 0.33203125 2.6308594 8.378906
11 1.3085938 1.8945312 2.480469 0.33203125 2.5714844 8.476562
12 1.6406250 1.8210937 4.980469 0.05859375 2.5925781 9.492188
13 2.7148438 3.0078125 4.023438 0.21484375 2.5753906 8.554688
14 0.7226562 1.7875000 3.085938 0.95703125 0.9734375 1.542969
Uses saved data files ripK_withcomments.csv for quadrat 47, I added columns for incidence and ripK_longR45.csv for quadrat 45, where I moved from a wide to a long format for quadrat 45 because it was easier to graph and then removed some values (NA) because at some time points there weren't really 2 peak clustering radii
ripK<-read.csv("ripK_withcomments.csv")
ripK$dates<-as.Date(ripK$dates,"%m/%d/%y")
ripK
dates p45.auc p47.auc bottom_upR45 upR45 top_upR45 bottom_lowR45
1 2019-03-04 12.79715782 5.1360509 3.6914062 4.316406 6.757812 0.8007812
2 2019-03-21 16.58548352 1.2491401 2.4804688 4.348438 10.000000 2.4804688
3 2019-04-11 0.01276313 1.1287722 7.3437500 7.422656 7.519531 NA
4 2019-05-02 10.76357877 9.5481176 0.7226562 5.320312 9.492188 NA
5 2019-05-16 4.04014681 7.3403452 NA NA NA 0.2539062
6 2019-05-28 3.07665302 11.0784857 3.6718750 4.397656 4.511719 0.2539062
7 2019-06-13 8.10033025 11.4838203 5.0781250 6.658594 10.000000 1.2304688
8 2019-07-01 7.05207728 9.7677461 5.1562500 5.912500 10.000000 1.2304688
9 2019-07-22 7.05207728 9.7677461 5.1562500 5.912500 10.000000 1.2304688
10 2019-08-16 3.35213051 9.8531471 6.1914062 6.855469 10.000000 1.6406250
11 2019-09-17 5.40229752 11.5345669 5.3320312 6.592969 10.000000 1.3085938
12 2019-10-14 4.02300920 15.1845803 5.2734375 6.625781 10.000000 NA
13 2019-11-12 3.98433812 5.9864486 5.3125000 5.882812 10.000000 2.7148438
14 2019-12-06 0.72366597 0.4679849 7.1679688 7.741406 9.082031 0.7226562
lowR45 top_lowR45 bottom_R47 R47 top_R47 p47.incidence p45.incidence
1 1.0625000 2.285156 0.80078125 1.5277344 2.929688 9 4
2 3.7296875 10.000000 1.30859375 1.7343750 2.753906 7 3
3 NA NA 2.44140625 3.8015625 6.132812 18 11
4 NA NA 1.01562500 4.5398438 10.000000 9 7
5 1.2609375 2.656250 0.33203125 3.2835938 10.000000 15 10
6 0.6078125 1.113281 0.60546875 3.6539062 9.375000 3 2
7 1.5859375 3.339844 0.33203125 3.1726563 8.906250 0 6
8 1.7031250 2.480469 0.33203125 3.0089844 9.042969 2 0
9 1.7031250 2.480469 0.33203125 3.0089844 9.042969 0 0
10 1.8945312 2.285156 0.33203125 2.6308594 8.378906 1 0
11 1.8945312 2.480469 0.33203125 2.5714844 8.476562 0 1
12 NA NA 0.05859375 2.5925781 9.492188 1 1
13 3.0078125 4.023438 0.21484375 2.5753906 8.554688 5 6
14 1.7875000 3.085938 0.95703125 0.9734375 1.542969 5 2
ripK$scaledI<-ripK$p47.incidence/10
ripK$p47.incidence
[1] 9 7 18 9 15 3 0 2 0 1 0 1 5 5
ripK$p45.incidence
[1] 4 3 11 7 10 2 6 0 0 0 1 1 6 2
topdifs<-ripK$top_R47-ripK$R47
bottomdifs<-ripK$R47-ripK$bottom_R47
ripK_p47<-ggplot(ripK)+
geom_point(mapping=aes(x=dates,y=R47,color=p47.auc))+
geom_path(mapping=aes(x=dates,y=R47),color="grey")+
geom_errorbar(aes(x=dates,ymin=R47-bottomdifs, ymax=R47+topdifs,color=p47.auc),width=3)+
scale_color_viridis(option = "magma",direction=-1,limits = c(0, 17))+
geom_bar(mapping=aes(x=dates,y=scaledI),stat="identity",fill=" dark grey")+
scale_y_continuous(name = "Cluster radius (m)", breaks=seq(0,10,1),expand = c(0, 0),limits=c(0,10),sec.axis = sec_axis(~.*10,name ="Disease incidence",breaks=c(0,4,8,12,16,20)))+
theme(legend.position = "none")+
theme(panel.background = element_blank())+
theme(text = element_text(family = "Times New Roman"),axis.text.x = element_text(angle = 70, hjust = 1))+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
labs(x="Survey Dates",color="ABC")+
scale_x_date( breaks=ripK$dates,date_labels="%d %b %y")+
guides(size = FALSE)
ripK_p47
longripK_45<-read.csv("ripK_longR45.csv")
longripK_45$dates<-as.Date(longripK_45$dates,"%m/%d/%y")
longripK_45
dates p45.auc bottom_R45 maxclust top_R45 measure p45.incidence
1 2019-03-04 12.79715782 3.6914062 4.3164062 6.757812 large 4
2 2019-03-21 16.58548352 2.4804688 4.3484375 10.000000 large 3
3 2019-04-11 0.01276313 7.3437500 7.4226563 7.519531 large 11
4 2019-05-02 10.76357877 0.7226562 5.3203125 9.492188 large 7
5 2019-05-16 4.04014681 NA NA NA large 10
6 2019-05-28 3.07665302 3.6718750 4.3976562 4.511719 large 2
7 2019-06-13 8.10033025 5.0781250 6.6585937 10.000000 large 6
8 2019-07-01 7.05207728 5.1562500 5.9125000 10.000000 large 0
9 2019-07-22 7.05207728 5.1562500 5.9125000 10.000000 large 0
10 2019-08-16 3.35213051 6.1914062 6.8554688 10.000000 large 0
11 2019-09-17 5.40229752 5.3320312 6.5929687 10.000000 large 1
12 2019-10-14 4.02300920 5.2734375 6.6257813 10.000000 large 1
13 2019-11-12 3.98433812 5.3125000 5.8828125 10.000000 large 6
14 2019-12-06 0.72366597 7.1679688 7.7414062 9.082031 large 2
15 2019-03-04 12.79715782 0.8007812 1.0625000 2.285156 small NA
16 2019-03-21 16.58548352 2.4804688 3.7296875 10.000000 small NA
17 2019-04-11 0.01276313 NA NA NA small NA
18 2019-05-02 10.76357877 NA NA NA small NA
19 2019-05-16 4.04014681 0.2539062 1.2609375 2.656250 small NA
20 2019-05-28 3.07665302 0.2539062 0.6078125 1.113281 small NA
21 2019-06-13 8.10033025 1.2304688 1.5859375 3.339844 small NA
22 2019-07-01 7.05207728 1.2304688 1.7031250 2.480469 small NA
23 2019-07-22 7.05207728 1.2304688 1.7031250 2.480469 small NA
24 2019-08-16 3.35213051 1.6406250 1.8945312 2.285156 small NA
25 2019-09-17 5.40229752 1.3085938 1.8945312 2.480469 small NA
26 2019-10-14 4.02300920 NA NA NA small NA
27 2019-11-12 3.98433812 2.7148438 3.0078125 4.023438 small NA
28 2019-12-06 0.72366597 0.7226562 1.7875000 3.085938 small NA
longripK_45$scaledI<-longripK_45$p45.incidence/10
longripK_45$p45.incidence
[1] 4 3 11 7 10 2 6 0 0 0 1 1 6 2 NA NA NA NA NA NA NA NA NA NA NA NA
[27] NA NA
ripK_p45<-ggplot(longripK_45)+
#geom_line(data = longripK_45[!is.na(longripK_45$maxclust), ],aes(x=dates,y=maxclust,group=measure),color="grey")+
geom_point(mapping=aes(x=dates,y=maxclust,color=p45.auc,group=measure,shape=measure))+
geom_bar(mapping=aes(x=dates,y=scaledI),stat="identity",fill="dark grey")+
geom_point(mapping=aes(x=dates,y=maxclust,color=p45.auc,group=measure,shape=measure))+
geom_path(mapping=aes(x=dates,y=maxclust,group=measure),color="grey")+
geom_errorbar(aes(x=dates,ymin=maxclust-(maxclust-bottom_R45), ymax=maxclust+(top_R45-maxclust),color=p45.auc),width=3)+
scale_color_viridis(option = "magma",direction=-1,limits = c(0, 17))+
scale_y_continuous(name = "Cluster radius (m)", breaks=seq(0,10,1),expand = c(0, 0),limits=c(0,10),sec.axis = sec_axis(~.*10,name ="Disease incidence",breaks=c(0,4,8,12,16,20)))+
theme(legend.position = "bottom")+
theme(panel.background = element_blank())+
theme(text = element_text(family = "Times New Roman"),axis.text.x = element_text(angle = 70, hjust = 1))+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
labs(x="Survey Dates",color="ABC",shape="Quadrat 45 R")+
scale_x_date( breaks=ripK$dates,date_labels="%d %b %y")+
guides(size = FALSE)
ripK_p45
ripK_p47/ripK_p45 +plot_annotation(tag_levels="a")
tiff("Figure7.tiff",width=180, height=220,units="mm",res=300)
ripK_p47/ripK_p45 +plot_annotation(tag_levels="A")
dev.off()
null device
1